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1  Introduction 


Whether  a  universal  quantum  computer  is  sufficiently  powerful  to  be  able  to  perform  quantum  field- 
theoretical  computations  efficiently  has  been  a  long-standing  and  important  open  question.  Efficient 
quantum  algorithms  for  simulating  quantum  many-body  systems  have  been  developed  theoretically 
m  and  implemented  experimentally  M,  but  quantum  field  theory  presents  additional  technical 
challenges,  such  as  the  formally  infinite  number  of  degrees  of  freedom  per  unit  volume.  In  earlier 
work  BIS],  we  presented  and  analyzed  a  quantum  algorithm  for  simulating  a  bosonic  quantum  held 
theory  called  </>4  theory.  That  algorithm  runs  in  a  time  that  is  polynomial  in  the  number  of  particles, 
their  energy,  and  the  desired  precision,  and  applies  at  both  weak  and  strong  coupling.  Hence,  it 
offers  exponential  speedup  over  existing  classical  methods  at  high  precision  or  strong  coupling.  In 
this  paper,  we  extend  our  work  to  fermionic  quantum  held  theories,  exemplihed  by  the  massive 
Gross-Neveu  model,  a  theory  in  two  spacetime  dimensions  with  quartic  interactions.  Although  our 
analysis  is  specihc  to  this  theory,  our  algorithm  can  be  adapted  to  other  massive  fermionic  quantum 
held  theories  with  only  minor  modification  while  retaining  polynomial  complexity. 

Our  quantum  algorithm  generates  scattering  events:  it  takes  (as  the  input)  the  momenta  of  the 
incoming  particles  and,  sampling  from  the  probability  distribution  of  possible  outcomes,  returns 
(as  the  output)  the  momenta  of  the  outgoing  particles  produced  by  the  physical  scattering  process. 
Physical  quantities  of  interest,  such  as  scattering  cross  sections,  can  thus  be  approximated  by 
repeated  runs  of  the  simulation,  together  with  statistical  data  analysis  similar  to  that  used  for 
particle-accelerator  experiments. 

The  features  of  fermionic  held  theories  not  present  in  bosonic  theories  pose  new  technical  prob¬ 
lems,  the  solutions  to  which  require  different  techniques.  Perhaps  the  most  obvious  difference  is 
the  anticommutation,  rather  than  commutation,  of  fermionic  helds.  This  forces  a  change  in  the 
representation  of  the  state  by  qubits:  we  use  an  encoding  method  for  fermionic  mode  occupation 
numbers  introduced  by  Bravyi  and  Kitaev  [9].  In  [8],  it  was  shown  that  simulation  of  Hamiltonian 
time  evolution  via  Suzuki- Trotter  formulae  has  efficiency  advantages  when  applied  to  spatially  local 
Hamiltonians.  Fermionic  anticommutation  makes  it  more  difficult  to  gain  efficiency  by  exploiting 
spatial  locality.  Nevertheless,  we  obtain  a  construction  that  gives  quasi-linear  asymptotic  scaling 
in  time  and  the  number  of  lattice  sites,  as  in  the  bosonic  case. 

In  contrast  with  bosonic  field  theories,  discretization  of  fermionic  field  theories  leads  to  the 
well-known  “fermion  doubling”  problem,  in  which  spurious  fermion  species  not  in  the  continuum 
theory  appear  in  the  discretized  theory.  One  solution  used  in  lattice  gauge  theory  is  to  add  to  the 
action  the  so-called  Wilson  term,  a  second-derivative  operator  that  vanishes  in  the  naive  continuum 
limit.  The  Wilson  term  can  also  be  accommodated  in  our  quantum  algorithm;  in  particular,  we 
show  how  it  can  be  turned  on  during  the  preparation  of  the  ground  state. 

In  general,  state  preparation  is  a  demanding  task.  The  algorithm  in  Bi  uses  a  three-step 
procedure.  First,  the  free  vacuum  is  prepared.  For  the  free  scalar  theory,  this  is  a  multivariate 
Gaussian  wavefunction.  Next,  wavepackets  are  excited  within  the  free  theory.  In  order  that  only 
single-particle  states  are  created,  an  ancillary  qubit  is  used,  together  with  a  particular  Hamiltonian 
that  acts  on  the  enlarged  space.  Finally,  the  interaction  is  turned  on  via  a  generalization  of  adiabatic 
state  preparation  that  can  be  applied  to  superpositions  of  eigenstates.  This  procedure  intersperses 
backwards  time  evolutions  governed  by  time-independent  Hamiltonians  into  the  turn-on  to  undo  the 
different  dynamical  phases,  which  otherwise  would  cause  undesirable  propagation  and  broadening 
of  wavepackets. 

The  state-preparation  method  analyzed  here  differs  from  that  of  urn  in  two  main  ways.  Prepa- 
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ration  of  the  free  vacuum  requires  modification  because  the  vacuum  of  the  free  fermionic  theory  is 
different  from  that  of  the  free  bosonic  theory.  For  this  purpose,  we  incorporate  a  separate  adia¬ 
batic  turn-on  step.  Furthermore,  sources  are  used  to  create  particle  excitations  after  the  coupling 
constant  is  adiabatically  turned  on,  rather  than  before.  (This  difference  is  not  required  by  the 
fermionic  nature  of  the  theory.)  This  method  has  the  advantage  that  it  works  when  bound  states 
are  possible,  in  which  case  the  adiabatic  wavepacket  preparation  of  might  fail.  Another  con¬ 
sequence  is  that  the  procedure  no  longer  requires  the  interleaving  of  backwards  time  evolutions  to 
undo  dynamical  phases.  On  the  other  hand,  a  disadvantage  is  that  the  preparation  of  each  particle 
has  a  significant  probability  of  producing  no  particle.  In  the  case  of  two-particle  scattering,  one  can 
perform  additional  repetitions  of  the  simulation,  and  recognize  and  discard  simulations  in  which 
fewer  than  two  particles  have  been  created.  However,  the  procedure  is  not  well  suited  to  processes 
involving  more  than  two  incoming  particles. 

We  analyze  two  different  measurement  procedures  to  be  used  as  the  last  step  of  the  simulation. 
The  first  method  is  to  return  adiabatically  to  the  free  theory  and  then  measure  the  number  operators 
of  the  momentum  modes.  For  unbound  states,  this  procedure  yields  complete  information  about 
particle  momenta,  but  is  not  well-suited  to  detecting  bound  states  or  resolving  spatial  information. 
The  second  procedure  is  to  measure  charge  within  local  regions  of  space.  These  measurements  can 
detect  charged  bound  states,  although  they  are  blind  to  neutral  ones.  Which  of  these  measurement 
schemes  is  preferable  depends  on  the  desired  application. 

There  is  a  substantial  body  of  work  on  analog  quantum  simulation  of  quantum  systems,  in¬ 
cluding  lattice  held  theories.  (See  (TO]  for  a  recent  review.)  In  such  work,  proposals  are  made  for 
the  engineering  of  experimental  systems  so  that  they  mimic  systems  of  interest,  that  is,  so  that 
the  Hamiltonians  of  the  laboratory  systems  approximate  Hamiltonians  of  interest.  The  proposed 
quantum  simulators  can  be  thought  of  as  specialized  quantum  computers.  In  contrast,  we  address 
digital  quantum  algorithms,  namely,  algorithms  to  be  run  on  a  universal,  fault-tolerant,  digital 
quantum  computer.  Our  work  thus  probes  the  fundamental  asymptotic  computational  complexity 
of  quantum  held  theories. 

There  is  also  an  extensive  literature  on  the  study  of  quantum  held  theories  on  classical  computers 
via  lattice  held  theory.  (See  Ch.  17  of  pT]  for  a  review  of  its  results  and  status.)  However,  classical 
lattice  algorithms  rely  on  analytic  continuation  to  imaginary  time,  t  — >•  —it.  Thus,  they  are  useful 
for  computing  static  quantities  such  as  mass  ratios,  but  are  unsuitable  for  calculating  dynamical 
quantities  such  as  scattering  cross  sections.  In  contrast,  our  quantum  algorithm  simulates  the 
dynamics  of  quantum  held  theories,  a  problem  that  is  expected  to  be  BQP-complete  and  thus 
impossible  to  solve  by  polynomial-time  classical  algorithms.  Although  our  algorithm  draws  upon 
some  concepts  from  lattice  held  theory,  new  techniques  are  needed,  particularly  for  state  preparation 
and  measurement. 

The  work  presented  in  this  paper  is  another  step  towards  the  goal  of  obtaining  an  efficient 
quantum  algorithm  for  simulating  the  Standard  Model  of  particle  physics.  Such  an  algorithm  would 
establish  that,  except  for  quantum-gravity  effects,  the  standard  quantum  circuit  model  suffices  to 
capture  completely  the  computational  power  of  our  universe. 

The  rest  of  this  paper  is  organized  as  follows.  Section  2  introduces  the  massive  Gross-Neveu 
model,  gives  an  overview  of  our  quantum  algorithm  for  computing  the  theory’s  scattering  ampli¬ 
tudes,  and  analyzes  the  algorithm’s  complexity.  Section  3  describes  in  detail  the  efficient  simulation 
of  the  Hamiltonian  time  evolution  in  the  quantum  circuit  model.  Section  4  presents  our  procedures 
for  state  preparation  and  measurement.  Finally,  Section  5  addresses  some  field-theoretical  aspects, 
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namely,  the  effects  of  a  non-zero  lattice  spacing  and  the  renormalization  of  mass,  which  are  crucial 
elements  in  our  complexity  analysis. 


2  Quantum  Algorithm 

In  this  section  we  describe  the  massive  Gross-Neveu  model  1H2.1I).  outline  the  steps  in  our  algorithm 
for  simulating  particle  scattering  processes  within  this  model  (H2.2I ),  and  give  an  overview  of  the 
algorithm’s  complexity  (H2.3I).  The  run  time  is  polynomial  in  the  inverse  of  the  desired  precision 
and  in  the  momenta  of  the  incoming  particles.  The  detailed  analysis  of  the  steps  of  the  algorithm 
that  contribute  to  the  overall  complexity  stated  in  H2.3I  is  given  in  later  sections. 


2.1  The  Massive  Gross-Neveu  Model 

The  theory  we  consider  is  a  generalization  of  the  Gross-Neveu  model  to  include  an  explicit  mass  term 
in  the  Lagrangian.  The  (original)  Gross-Neveu  model  [12]  is  a  quantum  field  theory  in  two  spacetime 
dimensions  consisting  of  N  fermion  species  with  quartic  interactions.  It  has  a  rich  phenomenology. 
Like  quantum  chromodynamics  (QCD),  the  theory  governing  the  strong  interactions,  it  has  the 
remarkable  property  of  asymptotic  freedom,  whereby  the  interaction  becomes  weaker  at  higher 
energies.  The  theory  has  a  discrete  chiral  symmetry,  ip  — >•  7 5ip,  where 
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This  symmetry  is  spontaneously  broken  by  the  non-perturbative  vacuum.  (The  related  theory 
known  as  the  chiral  Gross-Neveu  model  has  a  continuous  chiral  symmetry,  ip  — >•  e*®7  ip.)  Cor¬ 
respondingly,  mass  is  generated  dynamically,  and  the  theory  admits  a  topological  soliton,  the 
Callan-Coleman-Gross-Zee  (CCGZ)  kink.  Non-topological  solitons  also  exist  [13] . 

These  interesting  characteristics  have  attracted  intense  study  and  led  to  applications  not  only 
in  particle  physics  but  also  in  condensed-matter  physics,  including  studies  of  ferromagnetic  super¬ 
conductors  [14],  conducting  polymers,  and  systems  of  strongly  correlated  electrons  [15] . 

The  Gross-Neveu  model,  together  with  the  chiral  Gross-Neveu  model,  was  originally  solved 
in  the  limit  N  — >•  00  [12].  Via  inverse  scattering  methods  ns.  and  later  through  a  generalized 
Bethe  Ansatz  m  integrability  was  demonstrated  for  general  values  of  N,  a  feature  related  to  the 
existence  of  infinitely  many  conserved  currents  HH3-  The  model’s  5-matrix  is  factorizable  |T9ll20j: 
the  n-body  5-matrix  is  expressible  as  the  product  of  two-body  5-matrices. 

In  contrast,  the  massive  Gross-Neveu  model,  in  which  there  is  an  explicit  bare  mass,  is  thought 
not  to  be  integrable  for  arbitrary  values  of  N.  This  theory  still  exhibits  asymptotic  freedom,  but 
it  does  not  admit  solitons:  for  any  non-zero  mass,  the  CCGZ  kink  becomes  infinitely  massive  and 
disappears  [21] .  The  asymptotic  freedom  and  non-zero  bare  mass  make  a  rigorous  perturbative 
construction  of  the  theory  satisfying  the  Osterwalder-Schrader  axioms  possible  [22lf23] . 

The  massive  IV-component  Gross-Neveu  model  is  given  by  the  following  Lagrangian  in  two 
spacetime  dimensions: 
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where  each  field  ipj(x)  has  two  components,  7^  is  a  two-dimensional  representation  of  the  Dirac 
algebra,  and  ?/>  =  We  use  the  Majorana  representation,  namely, 
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The  components  of  the  field  operator  associated  with  the  particle  species  j  €  {1,  2, . . . ,  N}  will  be 
denoted  by  ipj}a,  ol  €  {0, 1}.  In  units  where  h  =  c  =  1,  any  quantity  has  units  of  some  power  of 
mass,  referred  to  as  the  mass  dimension.  We  shall  use  bold-face  to  represent  spatial  vectors,  such 
as  p  and  x,  to  distinguish  them  from  spacetime  vectors  x ^  =  ( t ,  x)  and  p ^  =  ( E ,  p).  Note,  however, 
that  we  are  considering  1+1  dimensions;  thus,  spatial  vectors  have  only  one  component. 

The  dimensionless  parameter  g  determines  the  strength  of  the  interaction.  When  g  =  0,  the  ijjj 
are  free  fields  obeying  the  Dirac  equation,  ( i'y^d L  —  mo)il>j(x )  =  0.  Then  one  can  write 


'ljjj(x) 


dp  1 

2tt  y/2 Kp 


(aj(p)u(p)e-iP-x  +  b](p)v(p)eiP-x)  , 


(4) 


where 

EP  =  \Jp2  +  ml ,  (5) 

Oj( p),  &j(p)  are  creation  and  annihilation  operators,  and  u,  v  satisfy 
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(m07°  -  P7VMP) 

=  -Epv{p) , 

(7) 

u\p)u(p)  =  ut(p)u(p) 

II 

to 

(8) 

u(p)t'u(-p) 

=  0, 

(9) 

h(p)u(p)  =  -v(pMp) 

=  2  m0  , 

(10) 

u(  pMp)  =  h(p)u(p) 

=  0. 

(11) 

In  the  Majorana  representation  (J3|),  one  has  the  following  concrete  solution: 
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2.2  Description  of  Algorithm 

To  represent  the  field  using  qubits,  we  first  discretize  the  quantum  field  theory,  putting  it  on 
a  spatial  lattice.  (Discretization  errors  are  analyzed  in  H5.ll)  Having  done  that,  our  algorithm 
consists  of  six  main  steps,  which  we  analyze  in  subsequent  sections. 

1.  Prepare  the  ground  state  of  the  Hamiltonian  with  both  the  interaction  term  (^q)  and  the 
nearest-neighbor  lattice-site  interactions  turned  off.  This  can  be  done  efficiently  because  the 
ground  state  is  a  tensor  product  of  the  ground  states  of  the  individual  lattice  sites. 

2.  Simulate,  via  Suzuki-Trotter  formulae,  the  adiabatic  turn-on  of  the  nearest-neighbor  lattice- 
site  interactions,  thereby  obtaining  the  ground  state  of  the  non-interacting  theory. 

1  The  Dirac  matrices  satisfy  {7+  7^}  =  7^7“'  +  7^7^  =  and  ipj{x)  is  a  spinor,  that  is,  its  Lorentz 

transformation  is  such  that  ©  is  Lorentz-invariant.  We  use  the  metric  (p’1'  =  diag(+l,  —1). 
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3.  Adiabatically  turn  on  the  interaction  term,  while  adjusting  the  parameter  mo  to  compensate 
for  the  renormalization  of  the  physical  mass. 

4.  Excite  particle  wavepackets,  by  introducing  a  source  term  in  the  Hamiltonian.  The  source 
term  is  chosen  to  be  sinusoidally  varying  in  time  and  space  so  as  to  select  the  desired  mass 
and  momentum  of  particle  excitations  by  resonance. 

5.  Evolve  in  time,  via  Suzuki- Trotter  formulae,  according  to  the  full  massive  Gross-Neveu  Hamil¬ 
tonian.  It  is  during  this  time  evolution  that  scattering  may  occur. 

6.  Either  use  phase  estimation  to  measure  local  charge  observables,  or  adiabatically  return  to 
the  free  theory  and  then  use  phase  estimation  to  measure  number  operators  of  momentum 
modes.  (The  choice  between  these  forms  of  measurement  depends  on  the  application.) 


2.3  Complexity 

In  this  section  we  bound  the  asymptotic  scaling  of  the  number  of  gates  needed  to  simulate  scattering 
processes  as  a  function  of  the  momentum  p  of  the  incoming  particles  and  the  precision  e  to  which 
the  final  results  are  desired.  The  effect  of  discretization,  via  a  lattice  of  spacing  a,  is  captured 
by  (infinitely  many)  terms  in  the  effective  Hamiltonian  that  are  not  present  in  the  continuum 
massive  Gross-Neveu  theory  fH5.ll).  Truncation  of  these  terms,  which  make  contributions  of  0(a) 
to  scattering  cross  sections,  therefore  constitutes  an  error.  Thus,  to  ensure  any  cross  section  a'  in 
the  discretized  quantum  field  theory  matches  the  continuum  value  a  to  within 

(1  —  e)a  <  a  <  (1  +  e)a,  (13) 


one  must  choose  the  scaling  a  ~  e  in  the  high-precision  limit,  that  is,  the  limit  e  — >■  0.  Similarly, 
in  the  large- momentum  limit,  one  must  choose  the  scaling  a  ~  p~4  in  order  to  ensure  that  the 
wavelength  of  each  particle  is  large  compared  with  the  lattice  spacing. 

It  suffices  to  use  an  adiabatic  process  of  duration 


T  =  0 
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(14) 


(where  L  is  the  length  of  the  spatial  dimension  and  m  is  the  physical  mass)  to  prepare  a  state 
within  a  distance  e  of  the  free  vacuum  (114.  II).  Using  Suzuki-Trotter  decompositions  of  the  form 
described  in  41.31  we  can  simulate  this  adiabatic  time  evolution  using  a  number  of  quantum  gates 
scaling  as 


G 


prep 


(15) 

(16) 


The  next  state-preparation  step  is  to  simulate  adiabatic  turn-on  of  the  coupling,  thereby  ob¬ 
taining  the  interacting  vacuum.  This  can  be  achieved  in  a  time  1H4.2I) 


Tturn— on 


o 


L2 


a4m3e 


(17) 
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Applying  Suzuki- Trotter  formulae,  one  obtains  a  gate  count  of 


G 


turn— on 


o 


L3 


a6m3e 


l+o(l)' 


(18) 


The  final  state-preparation  step  is  to  excite  particle  wavepackets.  We  do  this  by  applying  a 
time-dependent  perturbation  A W(t)  for  time  r.  It  is  necessary  to  choose  r  large  enough  and  A 
small  enough  to  suppress  the  production  of  particle  pairs.  The  choice  of  small  A  means  that  there 
will  be  a  substantial  probability  that  no  particle  is  produced.  Let  p\  denote  the  probability  that 
exactly  one  particle  is  produced.  In  a  typical  simulation  one  wishes  to  produce  an  initial  state  of 
two  spatially  separated  incoming  particles.  The  probability  that  both  of  these  are  produced  is  p\. 
The  simulations  in  which  one  or  both  initial  particles  has  failed  to  be  created  can  be  detected  at 
the  final  measurement  stage  of  the  simulation  and  discarded.  This  comes  at  the  cost  of  a  factor 
of  1/pi  more  repetitions  of  the  simulation.  The  probability  p\  is  independent  of  momentum  and 
scales  with  precision  as  p\  ~  e  144.31).  Also,  in  41. .11  one  finds  that  the  total  number  of  quantum 
gates  needed  for  the  excitation  step  is 

£— 4— o(i)  ,  as  e  -+  0  , 
p3+o(1)  ,  as  p  -+  oo  . 


In  both  the  high-momentum  and  high-precision  limits,  the  dominant  costs  in  the  algorithm  are 
the  two  adiabatic  state  preparation  steps,  whose  complexity  is  given  in  (1161)  and  (1181).  In  the  high- 
precision  limit,  to  compute  physical  quantities  such  as  scattering  cross  sections  to  within  a  factor  of 
(1  +  e),  one  must  choose  a  to  scale  as  e  1 115. II).  Also,  in  this  limit,  the  complexity  contains  a  further 
factor  of  1  /e  owing  to  postselection  of  simulations  in  which  both  wavepacket  excitations  have  been 
successful  144.31).  Substituting  a  ~  e  into  (fT6j)  and  including  this  extra  factor  of  1/e  yield  a  total 
complexity  of  ©(e-8-0^1').  In  the  high-momentum  limit,  a  must  scale  as  1/p  to  ensure  that  the 
particle  wavelength  is  long  compared  to  the  lattice  spacing,  and  L  must  scale  as  p  to  accommodate 
the  excitation  step  144.3ft .  In  summary,  we  obtain 


(  0(e  8  °(1)) ,  as  e  ->  0 , 

\  O(p9+°(0) ,  as  p  +  oo  . 


(20) 


Note  that  these  are  only  upper  bounds  on  the  complexity,  and  it  may  be  possible  to  improve  them 
by  using  more  detailed  analysis,  such  as  more  specialized  adiabatic  theorems. 


3  Qubits  and  Quantum  Gates 

We  divide  the  problem  of  simulating  Hamiltonian  time  evolutions  in  the  massive  Gross-Neveu  model 
into  three  subproblems.  The  first  subproblem  is  to  represent  the  state  of  the  field  with  qubits.  We 
do  this  by  choosing  a  complete  set  of  commuting  observables  and  encoding  their  eigenvalues  with 
strings  of  bits  143.11).  The  second  subproblem  is  to  simulate  local  fermionic  gates  on  the  degrees  of 
freedom  defined  by  the  commuting  observables.  Achieving  this  in  an  efficient  manner  is  non-trivial 
because  of  the  fermionic  statistics.  For  this  purpose,  we  employ  a  technique  due  to  Bravyi  and 
Kitaev  [9],  which  implements  fermionic  statistics  with  only  logarithmic  overhead  in  the  number 
of  lattice  sites  143.21).  The  third  subproblem  is  to  decompose  the  time  evolution  governed  by  the 
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massive  Gross-Neveu  Hamiltonian  into  a  product  of  local  fermionic  gates.  We  do  this  using  high- 
order  Suzuki- Trotter  formulae  [53]  with  optimizations  tailored  to  the  fermionic  statistics  and  the 
spatially  local  nature  of  the  Hamiltonian  HI3.3I).  The  local  unitary  transformations  act  on  at  most 
2 2 ,v - d i m e n s i o n a  1  Hilbert  spaces  and  can  therefore  be  efficiently  decomposed  into  elementary  gates 
for  any  constant  number  of  particle  species,  N,  via  the  Solovay-Kitaev  algorithm  [251126] . 

3.1  Representation  by  Qubits 

First,  we  put  the  massive  Gross-Neveu  model  on  a  spatial  lattice 

n  =  aZL.  (21) 

For  simplicity,  we  impose  periodic  boundary  conditions,  so  that  H  can  be  considered  a  circle  of 
circumference  L  =  aL.  The  Hamiltonian  is 


H  =  H0  +  Hq  +  Hw  , 


where 


N 


Hq  = 


xEQ  j= 1 


.  l  ipj  (x  +  a)  —  if}j  (x  —  a)  . 

-*7  — - 2^ - hmoV’j(x) 


H„  =  - 


9o 


N 


Hw  = 


2  5Za( 

xEQ  ' j= 1 
N 

[_^^x)  Wx  +  a)  ~  2^i(x)  +  Vb(x  -  a)) 


xgf2  j=l 


(22) 


(23) 

(24) 

(25) 


Here,  Hg  is  the  interaction  term,  and  Hw  is  the  Wilson  term,  used  to  prevent  fermion  doubling 
m-  Correspondingly,  0  <  r  <  1  is  called  the  Wilson  parameter.  H  is  spatially  local  in  the  sense 
that  it  consists  only  of  single-site  and  nearest-neighbor  terms  on  the  lattice. 

Let  r  denote  the  momentum-space  lattice  corresponding  to  H,  namely, 

(26) 

We  can  deduce  the  spectrum  Hq  +  Hw  using 

Vt(x)  =  ^(°i(pHp)e<p'x  +  65(p)v(p)e“<p'x)  ,  (27) 

per  L  V2Ep  V  7 

&(x)  =  E7i(fl!(P)“(P)e'iP'X  +  ^(P)“(P)ei'  (28) 

per  L  V2Ep  V  7 

The  inverse  transformation  is 

aj (p)  =  SE  °e-»p-x^-(x) ,  (29) 

V  xeo 

6](p)  =  (30) 

V2EP  xen 
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Substituting  ([271)  and  (f28l)  into  (1231)  and  (1251)  and  neglecting  the  vacuum  energy,  we  obtain 


N 

Ho  +  Hw  =  EEi4“>  (m0)  (a\{p)aj(p)  +  &](p)6j(p))  , 
3= i per 

where 

Ep\m0)  =  \J  (m0  +  ^  sin2  +  ^  sin2  (pa). 

From  the  canonical  fermionic  anticommutation  relations 

=  Wj>(X)>^fe,^(y)}  =  °: 

it  follows  that 


{oj(p),4(q)}  =  LSPtqSjjkt, 

{bj(p),bl(q)}  =  LSP}Cl5jjkl, 


(31) 


(32) 


(33) 

(34) 

(35) 

(36) 


with  all  other  anticommutators  involving  a  and  b  operators  equal  to  zero.  We  thus  have  the 
following  interpretation:  there  are  N  independent  fermion  species,  created  (with  momentum  p) 
by  a|(p), . . .  ,ajy(p)  and  annihilated  by  ai(p), . . . ,  ajv(p)-  Similarly,  for  each  species  j,  bj( p)  and 
bj( p)  are  the  creation  and  annihilation  operators  for  a  corresponding  antifermion.  Thus,  H  acts 
on  a  Hilbert  space  of  dimension  22NL . 


We  can  specify  a  basis  for  the  Hilbert  space  of  field  states  by  choosing  a  complete  set  of 
commuting  observables.  The  basis  is  then  indexed  by  the  set  of  eigenvalues  of  these  observables. 
The  fermionic  anticommutation  relations  {a,  a^}  =  1,  {a,  a}  =  0  imply  that  the  algebra  generated 


by  a  and  a t  has  the  irreducible  representation  a  — » 


0  1 
0  0 


0  0 
1  0 


which  is  unique  up 


to  the  choice  of  basis.  Hence,  the  eigenvalues  of  a^a  are  0  and  1.  The  two  basis  vectors  for  the 
space  on  which  a  and  a)  act  are  interpreted  as  the  presence  or  absence  of  a  fermion. 

Thus,  by  (f33l)  and  (1341). 


S3 


=  {aifi 


(x)V’j.a(x) \j  =  1,...,N;  a  =  0, 1;  x  G  H} 


(37) 


is  a  set  of  2NL  commuting  observables,  each  of  which  has  eigenvalues  zero  and  one.  Similarly,  by 
(1351)  and  (f36l). 

Sp  =  {L~1aij(p)aj(p)\j  =  1,...,N;  p  €  T}  U  {L~1bj(p)bj  (p)\j  =  1, . . .  ,  TV;  p  €  T}  (38) 


is  a  set  of  2NL  commuting  observables,  each  with  eigenvalues  zero  and  one.  In  the  non-interacting 
theory,  the  eigenvalues  of  the  elements  of  Sp  are  interpreted  as  the  fermionic  occupation  numbers 
of  different  momentum  modes. 

The  Hamiltonian  Ho  +  H\y  is  called  the  free  theory.  The  eigenstates  of  the  number  operators 
in  Sp  are  eigenstates  of  Ho  +  H\y ,  and  thus  the  particles  do  not  interact.  The  rest  mass  of  these 
non-interacting  particles  is  Eq  (mo)  =  mo-  It  is  not  known  how  to  solve  for  the  spectrum  of 
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Hq  +  H\v  +  Hg  analytically,  but  the  eigenvalue  spectrum  of  Hq  +  H\y  +  Hg  can  still  be  characterized 
in  terms  of  particles.  The  rest  mass  m  of  the  particles  in  Ho  +  H\y  +  Hg  is  equal  to  the  eigenvalue 
gap  between  the  ground  state  (also  called  the  vacuum)  and  the  first  excited  state.  In  the  interacting 
theory,  it  is  no  longer  true  that  m  =  mo-  Rather,  m  depends  in  a  non-trivial  way  on  mo,  go,  and 
a;  the  mass  is  said  to  be  renormalized.  A  quantitative  analysis  of  this  effect  contributes  to  our 
analysis  of  adiabatic  state  preparation  and  is  given  in  H5.2I 

One  can  represent  the  quantum  state  of  the  fermionic  fields  using  2 NL  qubits  to  store  the 
eigenvalues  of  the  elements  of  either  Sx  or  Sp.  The  ground  state  of  the  free  theory  in  the  Sp 
representation  is  thus  1 000 . . .).  However,  the  ground  state  of  the  interacting  theory  is  non-trivial 
in  both  representations.  We  define  our  qubit  basis  in  terms  of  the  elements  of  Sx,  because  the 
Gross-Neveu  Hamiltonian  is  local  in  this  basis,  which  improves  the  scaling  of  the  Suzuki-Trotter 
formulae  used  to  implement  time  evolution.  However,  we  do  not  simply  store  the  eigenvalues  of 
the  elements  of  Sx  directly  as  the  values  of  the  qubits.  This  representation  would  be  somewhat 
inefficient  to  act  upon,  because  direct  implementation  of  the  fermionic  minus  signs  requires  0(L ) 
gates.  Instead,  we  apply  the  method  of  [9]  to  reduce  this  overhead  to  O(logT),  as  described  next. 

3.2  Simulating  Fermionic  Gates 

The  implementation  of  fermionic  gates  using  qubits  can  present  a  technical  challenge  [9].  As  an 
example,  consider  the  unitary  transformation  Uj:0l(x.)  =  + -i/>ja(x)).  This  toggles  the 

eigenvalue  of  aVh,a(x)V(}  Q(x)  between  zero  and  one.  Such  a  toggling  can  be  implemented  on  qubits 
with  the  NOT  gate.  However,  to  satisfy  the  fermionic  anticommutation  relations  (13311  and  (13411  the 
sign  of  the  transition  amplitude  between  the  zero  and  one  state  must  depend  on  the  occupation  of 
other  modes.  A  well-known  way  to  satisfy  ()33(i  and  (I34|)  is  to  use  a  Jordan-Wigner  transformation, 
in  which  the  modes  are  given  an  ordering  and  C/jiQ,(x)  is  represented  by  the  operator  ax®az®. .  .®az, 
where  the  az  operators  apply  to  all  preceding  mode^l  [28].  Unfortunately,  this  method  clearly  has 
an  0(L )  overhead.  In  [9],  Bravyi  and  Kitaev  give  a  method  with  only  O(logL)  overhead,  which  we 
briefly  review  here. 

Let  rii  be  the  occupation  number  of  the  zth  fermionic  mode  according  to  some  chosen  numbering 
of  the  modes  from  1  to  2NL.  To  implement  the  minus  signs  in  UjiQ(x),  one  needs  to  know 
where  the  sum  is  over  all  preceding  modes.  Thus,  a  natural  encoding  of  fermionic  mode  occupation 
numbers  is  to  store  the  quantities  L  =  X)’=i  rij  instead  of  the  quantities  rij.  This  encoding  has  the 
advantage  that  calculating  the  relevant  signs  has  an  0(1)  cost.  However,  it  has  the  disadvantage 
that,  if  the  occupation  number  of  the  ith  mode  changes,  then  i  —  1  of  the  U  values  must  be  updated. 
Thus,  updates  have  an  O(L)  cost.  The  Bravyi-Kitaev  encoding  uses  the  following  compromise,  in 
which  the  calculation  of  the  relevant  signs  and  the  update  steps  can  both  be  performed  in  time 
O(logL). 

The  mode  index  i  €  {1, . . . ,  2 NL}  can  be  represented  by  a  bit  string  of  length  l  =  (log2(2ArL)] . 
One  can  define  the  following  partial  order  on  these  bit  strings.  Consider  two  bit  strings  x  = 
xixi^i  . . .  x\  and  y  =  yiyi-i  ■  ■  ■  y\.  Then  x  A  y  if,  for  some  r,  Xj  =  yj  for  j  >  r  and  yr-\  =  yr- 2  = 
. . .  =  yi  =  1.  Now,  let  kj  =  n«-  ^ny  total  occupation  number  tt  can  be  computed  from  the 

kj  quantities  in  0(log  L )  time  and  changing  the  occcupation  of  any  mode  rij  requires  updating  only 
O(logL)  of  the  kj  quantities  [9]. 

2Note  that  one  can  apply  both  the  Jordan-Wigner  and  Bravyi-Kitaev  methods  for  implementing  fermionic  oper¬ 
ators  on  quantum  computers  in  any  number  of  spatial  dimensions,  using  an  arbitrary  numbering  of  modes. 
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In  fact,  the  Bravyi-Kitaev  construction  is  relevant  only  to  the  excitation  of  wavepackets  1114.311. 
In  all  other  parts  of  our  algorithm,  we  simulate  a  Hamiltonian  in  which  every  term  is  a  product  of 
an  even  number  of  fermionic  field  operators,  all  acting  on  the  same  site  or  on  nearest-neighbor  sites 
in  one  dimension.  In  this  case,  traditional  Jordan-Wigner  techniques  incur  only  0(1)  overhead. 


3.3  Application  of  Suzuki- Trotter  Formulae  to  Fermionic  systems 

In  this  section,  we  describe  how  to  construct  efficient  quantum  circuits  that  simulate  time  evolution 
induced  by  the  Hamiltonian  H  defined  in  ([22]).  ([23]) .  ([24]) .  and  ([25]) .  We  present  the  case  in  which 
H  is  time- independent.  By  the  results  of  [29],  the  same  analysis  applies  to  the  simulation  of  the 
time-dependent  Hamiltonians  that  we  use  in  adiabatic  state  preparation.  (See  also  [30].) 

Using  a  fcth-order  Suzuki- Trotter  formula,  one  can  implement  Hamiltonian  time  evolution  of 
duration  t  using  a  number  of  quantum  gates  that  scales  as  f1+ ^  [24[[3T] .  Generally,  applying  a 
Suzuki- Trotter  formula  directly  to  a  Hamiltonian  of  the  form 

m 

H  =  YJHi  (39) 

i—1 

yields  an  algorithm  with  0(m1+o^)  timesteps,  and  hence  0(m  s+oM)  gates,  if  the  Hi  are  not 
mutually  commuting.  Thus,  it  is  often  advantageous  to  group  terms  in  a  Hamiltonian  like  (1391)  into 
as  small  a  collection  as  possible  of  sets  of  mutually  commuting  terms  [11132]. 

Consider  the  problem  of  simulating  the  Hamiltonian  H  defined  in  (|22l) ,  (1231) ,  ()24D ,  and  (1251) .  By 
(1331)  and  (IMl).  one  sees  that 

fe(x)^-(x),V'fc(y)V,fc(y)]  =o,  (40) 

regardless  of  whether  j  =  k  or  x  =  y.  Thus,  we  start  by  decomposing  H  as  a  sum  of  two  parts, 
the  single-site  terms  and  the  terms  that  couple  nearest  neighbors: 


H  =  HSS  +  Hn 


(41) 


where 


H ss  —  ^  '  o, 

xgO  Lj=l 


N  2  /  N 

Y  (mo^(x)^(x)  +  ^(x)^-(x))  +  y  (  Xl^^x)V’j(x) 
3  =  1  v  3= 1 


By  (1401).  e  iHssSt  decomposes  into  a  product  of  local  unitary  transformations. 
All  terms  in  Hnn  are  of  the  form 


(42) 


^la(X)Vh,/3(y)  +  4/3(y)^>(x)  .  (43) 

for  x  =  y  ±  a.  Terms  with  a  =  j3  and  terms  with  are  both  present  in  Hnn. 

Given  an  operator  of  the  form  (1431).  let  us  refer  to  the  subset  of  {1, . . . ,  N}  x  {0, 1}  x  SI  on  which 
it  acts  as  its  support.  Because  they  consist  of  a  product  of  an  even  number  of  fermionic  operators, 
any  two  operators  of  the  form  ()43l)  commute  provided  they  have  disjoint  support.  Thus,  we  next 
decompose  Him  as 

Hnn  =  Hi  +  H2  +  H3  +  Hi  ,  (44) 


where  each  of  H\, ... ,  H±  consists  of  a  sum  of  terms  with  non- intersecting  support. 
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Figure  1:  Vertices  represent  elements  of  {0, 1}  x  Q  two  vertices  are  connected  by  an  edge  if  Hnn  couples  these 
sites.  (Different  species  are  never  coupled  by  Hn n,  so  the  full  graph  with  vertices  corresponding  to  elements 
of  {1, ... ,  N}  x  {0,  ljxll  would  consist  of  N  disconnected  copies  of  the  graph  shown.)  The  edges  can  be 
colored  with  four  colors  such  that  each  node  has  no  more  than  one  incident  edge  of  each  color.  One  can 
obtain  the  decomposition  Hnn  =  Hi  +  Hi  +  +  H4  by  choosing  H 1  to  be  the  sum  of  all  interaction  terms 

along  the  edges  labeled  1  (which  are  blue),  H2  to  be  the  sum  of  all  the  interaction  terms  along  edges  labeled 
2  (which  are  red),  and  so  on. 


In  Hnn  there  is  no  coupling  between  different  species,  that  is,  no  products  of  ipj  and  for 
j  ^  k.  Thus,  we  can  ignore  the  index  j.  We  now  construct  a  graph  whose  vertices  correspond  to 
the  elements  of  {0, 1}  x  12.  We  draw  an  edge  between  two  vertices  if  there  exists  a  term  in  Hnn 
with  the  corresponding  support.  One  sees  that  this  graph  is  as  shown  in  Fig.  |T]  The  graph  is 
edge-colorable  with  four  colors,  and  therefore  Hnn  is  correspondingly  decomposable  as  in  (|44|)  with 
each  of  H\ .  H2,  H%,  Hi  consisting  of  a  sum  of  commuting  terms.  (Because  of  the  periodic  boundary 
conditions,  this  works  only  if  L  is  even,  which  we  assume  henceforth.) 

The  unitary  time  evolution  induced  by  H  =  Hss  +  H\  +  H2  +  H$  +  Hi  can  be  approximately 
decomposed  via  high-order  Suzuki- Trotter  formulae  into  a  sequence  of 

ns_T  =  0((t/a)1+o(1)Lo(1)e_o(1))  (45) 

time  evolutions  induced  by  individual  members  of  {Hss,  Hi,  K/2,  #3,  H. 4}.  The  scaling  with  t  follows 
from  [241129] .  The  scaling  with  L  is  a  consequence  of  the  spatial  locality  of  H  (see  §4.3  of  [£]),  that 
is,  the  property  that  only  nearest-neighbor  sites  are  coupled.  The  scaling  with  a  is  a  consequence 
of  the  fact  that  the  individual  terms  in  the  Hamiltonian  each  have  norm  at  most  of  order  a-1.  This 
affects  the  magnitude  of  the  error  term  in  the  Suzuki- Trotter  decomposition,  which  arises  from 
commutators  of  these  terms. 

Each  member  of  {Hss,  H 1,  H2,  H3,  H4}  is  a  sum  of  0(L )  commuting  terms.  The  time  evolution 
induced  by  commuting  terms  Mi  decomposes  as  If  each  Hi  acts 

on  only  a  constant  number  of  qubits,  then  the  individual  factors  e~lHit  in  this  product  can  each 
be  simulated  in  0(1)  time,  by  the  Solovay-Kitaev  theorem  [251(26],  Thus,  including  a  logarithmic 
overhead  for  fermionic  statistics,  the  cost  of  implementing  e~lJt  for  any  J  £  {Hss,  Hi,  H2,  H3,  H4} 
is  O(L).  By  ([45]).  the  total  cost  of  time  evolution  is  0(  (^)1+°^  e“°A))  quantum  gates. 

4  State  Preparation  and  Measurement 

We  divide  the  problem  of  state  preparation  into  three  steps,  described  in  44. 11-44.31  preparing  the 
free  vacuum,  transforming  the  free  vacuum  into  the  interacting  vacuum,  and  exciting  wavepackets 
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on  the  background  of  the  interacting  vacuum.  Two  possible  measurement  procedures  are  described 
in  41.11  and  84.51 

4.1  Preparing  the  Free  Vacuum 

Although  the  free  Hamiltonian  H$  +  Hw  is  exactly  solvable,  preparing  its  ground  state  in  the  Sx 
representation  on  a  quantum  computer  is  non-trivial.  We  do  so  using  adiabatic  state  preparation, 
as  follows.  Let 


H(s)  = 


E 


N 

aX^'( 

3= 1 


X 


-szq 


1 4’j  (x  +  a)  -  ipj  (x  -  a) 


2  a 


+  mifjfx) 


+  sH\y. 


(46) 


The  energy  gap  of  this  Hamiltonian  is  equal  to  the  parameter  m  for  all  s.  We  set  this  equal  to  the 
physical  mass  of  the  particles  whose  scattering  we  ultimately  wish  to  simulate. 

14(0)  is  a  sum  of  separate  Hamiltonians  acting  on  each  lattice  site  and  each  species  of  parti¬ 
cle.  Its  ground  state  is  therefore  the  tensor  product  of  the  ground  states  of  the  four-dimensional 
Hilbert  spaces  associated  with  each  pair  (x, j)  £  x  {1,...,1V}.  (Specifically,  the  ground  state 
for  a  given  site  is  -^=  (|01)  +  i|10}),  where  \bobi)  with  bo , b i  €  {0,1}  denotes  the  state  satisfying 

a,i/’jo(x)'0j,o(x)|6o&i)  =  bo\bobi)  and  aV’j  1(x)^ii(x)|6o&i)  =  bfibobi).)  The  cost  of  producing  this 
tensor  product  of  NL  local  states,  including  the  cost  of  fermionic  antisymmetrization  via  the  en¬ 
coding  of  [9],  is  0{NL\og{NL)). 

After  the  ground  state  of  Hq  has  been  prepared,  the  complexity  of  the  remaining  adiabatic  state 
preparation  is  determined  by  the  adiabatic  theorem  [MIEU . 


Theorem  1.  Let  H{s)  be  a  finite- dimensional  twice  differentiable  Hamiltonian  on  0  <  s  <  1  with 
a  non- degenerate  ground  state  |0o(s))  separated  by  an  energy  gap  'y(s).  Let  | ip(t)}  be  the  state 
obtained  by  Schrodinger  time  evolution  according  to  the  Hamiltonian  H(t/T)  from  the  state  |^>o(0)) 
at  t  =  0.  Then,  with  an  appropriate  choice  of  phase  for  \(f>o (t)),  the  error  A  =  ||  \  fi(T))  —  |</>o(l))  || 
satisfies 


A 


1 

lW 


d  H 
ds 


1 


d  H 
ds 


d  H 

2  1 

H - rr 

d214 

)1 

/o  y73 

ds 

72 

ds2 

>\ 

(47) 


Analyzing  the  adiabaticity  of  this  process  is  relatively  easy,  because  (l27|)  and  (1281)  diagonalize 
H(s)  (and  ^)  for  all  s.  One  finds  that  the  eigenvalue  gap  of  H(s)  throughout  the  adiabatic  path 
0  <  s  <  1  is  always  precisely  m.  Furthermore, 


dH 

ds 


fiEp] ’(°)  (aj(p)ai(p) 

l=i  per 


+  &3(p)Mp))  • 


(48) 


Thus, 


dH 

ds 


2N^2E{pa)(0) . 

per 


(49) 
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For  large  L,  ^)pgr  X  becomes  well  approximated  by  the  integral  JQ27r,/adp.  Thus,  using  (1321).  we 
obtain 


dH 

ds 


(50) 

(51) 

(52) 


where 


r;(r)  =  J  dp]J 4 r2  sin4  +  sin2  ( p ) . 


(53) 


We  can  therefore  substitute  =  0,  ||^||  =  0(La  2)  and  7  =  m  into  (ITH).  Theorem  [T]  then 
shows  that  we  can  prepare  a  state  with  distance  no  more  than  eprep  from  the  exact  state  using 


T  =  0 


L2 


a^m^e 


prep 


(54) 


Note  that  the  adiabatic  theorem  applied  here,  though  convenient  because  of  its  generality,  may  not 
yield  a  tight  upper  bound  on  the  run  time. 


4.2  Preparing  the  Interacting  Vacuum 

Given  the  ground  state  of  the  free  theory,  we  can  prepare  the  ground  state  of  the  interacting 
theory  by  adiabatically  varying  the  parameters  g q  and  mo  in  the  massive  Gross-Neveu  Hamiltonian, 
starting  from  g q  =  0.  For  adiabaticity  to  be  maintained,  the  physical  mass  must  not  vanish  at  any 
point  in  the  adiabatic  path.  By  45.21  the  physical  mass  varies  with  g p  according  to 


m  =  m0-  cigl  -  c2go  +  0(gl) , 

(55) 

where  c\ ,  C2  >0  are  given  by 

Cl 

m  ,  /  1  \ 

=  log  ( - )+•••, 

Z7T  \ma / 

(56) 

C2 

~  ra  (9.3V  —  0.07)  log2(mo)  +  •  •  •  . 

l07 T6 

(57) 

(The  coefficients  in  (157jl  were  determined  numerically.)  The  vanishing  of  the  physical  mass  marks 
the  location  of  a  quantum  phase  transition,  which  cannot  be  adiabatically  crossed.  Equation  (15511 
indicates  that  the  phase  diagram  takes  the  schematic  form  as  shown  in  Fig.  [2j 

As  in  <14.11  we  parametrize  our  adiabatic  state  preparation  by  a  quantity  s ,  which  increases  over 
time  from  0  to  1.  In  this  second  adiabatic  process,  the  Hamiltonian  is  the  full  massive  Gross-Neveu 
Hamiltonian  with  s-dependent  parameters  #o(s)  and  mo(s).  We  choose  5g(0)  =  0  and  mo(0)  =  m  so 
that  the  initial  Hamiltonian  of  this  adiabatic  process  matches  the  final  Hamiltonian  of  the  preceding 
adiabatic  step.  Thus,  the  ground  state  at  s  =  0  is  the  free  vacuum  prepared  in  the  previous  step 
of  the  algorithm.  To  keep  our  analysis  simple,  we  choose  a  linear  adiabatic  path,  namely, 

ffo(«)  =  s9q  , 

m0(s)  =  m  +  s6m.  (58) 
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Figure  2:  Our  perturbative  calculations  of  the  physical  mass  in  the  massive  Gross-Neveu  model  indicate  a 
phase  diagram  with  the  qualitative  features  illustrated  above.  The  phase  above  the  dashed  curve  is  accessible 
adiabatically  from  the  free  theory  but  the  phase  below  is  not.  The  arrow  depicts  our  linear  adiabatic  path, 
described  in  (15911 .  Our  perturbative  analysis  shows  that  the  first  two  derivatives  of  the  phase  transition  curve 
with  respect  to  are  both  positive  and  diverge  only  as  poly(log(moa))  in  the  limit  a  — >  0. 


We  choose  5m  so  that  the  physical  mass  at  s  =  1  is  equal  to  the  physical  mass  at  s 
order  in  g q, 

=  ci9o  +  c2fifo  +  '  ”  ) 


as  illustrated  in  Fig.  [2j 
By  (PI),  =  0  and 


0.  To  second 
(59) 


dH 

ds 


£ 

xSO 


9o 


N 

X^7x)Vh(x) 

i=i 


(60) 


Furthermore,  the  minimal  eigenvalue  gaps  occur  at  s  =  0  and  s  =  1  and  are  equal  to  the  final 
physical  mass  m.  Thus,  to  apply  Theorem Q]  we  need  only  bound  ||^||. 

We  can  deduce  the  spectrum  of  ^  by  the  following  transformation: 


«i(x) 


1 

71 

i 

71 


{ipj, o(x)  -  iipj, i(x)) 

(^,o(x)  +  #j,l(x)) 


This  corresponds  to 

Vh(x)  =  y=  (aj(x)u(O)  +  &5(x)w(0))  , 
where  u,  v  are  defined  in  (fl2l).  Using  (1331)  and  (l34l).  one  can  verify  that 


(61) 

(62) 

(63) 


(ai(x),4(y)}  =  {^(x),4(y)}  =  a  Hj^kSx, yl,  (64) 

(ai(x),afe(y)}  =  (6j(x),6fe(y)}  =  0,  (65) 

(aj(x)A(y)}  =  (at(x),6fc(y)}  =  0.  (66) 


Thus,  (ij  (x) .  (x) ,  bj  (x) ,  b-  (x)  are  creation  and  annihilation  operators  for  2 N  species  of  fermions 
localized  on  the  spatial  lattice.  By  (|63l). 


^■(x)^'(x)  =  7(x)ai(x)  _  Mx)kj(x)  > 


(67) 
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from  which  we  obtain 


N 


X^(x)^'(x) 

3= 1 


=  2Na 


-i 


and  hence 


dH 


ds 


2Lg2N2 


—  5m2NL  + 


By  the  results  of  135.21  we  find  that  5m  =  0(log2(ma)).  Hence,  recalling  that  L 
obtain 


dH 


ds 


Therefore,  by  Theorem  []]  the  diabatic  error  is  at  most 

e  =  O 

=  O 


1  II  dH  | 

^  II  ds  I 


^turn-on  7^ 

L2 


7"turn— 


It  thus  suffices  to  choose 


Tfiirn — on  —  O 


L 2 


a^em3 


(68) 

(69) 
L/a ,  we 

(70) 

(71) 

(72) 

(73) 


In  the  above  procedure,  we  choose  our  adiabatic  path  so  that  the  initial  and  final  physical  masses 
equal  some  user-specified  value  m.  To  achieve  this,  one  needs  to  tune  the  quantity  5m  in  accordance 
with  1158])  and  (l59jl.  For  sufficiently  weak  coupling,  the  proper  choice  of  5m  can  be  determined  by 
the  perturbative  calculations  performed  in  115.21  In  the  strongly  coupled  case,  these  perturbative 
calculations  no  longer  provide  precise  guidance  as  to  a  choice  of  5m.  Instead,  as  previously  discussed 
in  [8],  the  adiabatic  path  can  be  determined  by  the  quantum  computer.  Specifically,  one  can 
measure  the  physical  mass  at  a  given  coupling  strength  go  by  exciting  a  particle  and  measuring 
energy  via  phase  estimation.  This  measurement  guides  the  choice  of  a  suitable  adiabatic  path  to 
a  slightly  larger  coupling  strength,  at  which  point  the  mass  can  be  measured  again.  Iterating  this 
process,  one  can  reach  any  coupling  strength  for  which  the  corresponding  vacuum  is  in  the  same 
quantum  phase  as  the  free  vacuum. 


4.3  Exciting  Wavepackets 

After  preparing  the  interacting  vacuum,  |vac),  we  excite  wavepackets  by  simulating  a  source  that 
varies  sinusoidally  in  space  and  time  so  as  to  induce  excitations  of  some  particular  total  energy  and 
momentum  by  resonance.  Given  the  physical  rest  mass  m  of  the  particles,  we  can  choose  this  energy 
and  momentum  so  that  the  only  corresponding  state  is  a  single-particle  state.  (For  a  given  total 
momentum,  an  unbound  state  of  two  particles  will  have  greater  energy  than  the  corresponding  state 
of  one  particle.  In  the  ultrarelativistic  limit,  p>m,  this  energy  difference  scales  as  m2/p.)  In  the 
remainder  of  this  section,  we  show  that,  using  a  source  of  spatial  extent  l  and  duration  r,  one  can 
ensure  that  excitations  off  resonance  are  suppressed  as  ~  exp  [—  |  (Z2( p  —  po)2  +  t2(E  —  Eb)2)] . 
Hence,  by  simulating  a  process  of  duration  r  ~  p/m 2  and  spatial  extent  l  ~  p/m2,  one  can  control 
the  incoming  momentum  and  ensure  that  the  probability  of  obtaining  more  than  one  particle  is 
small. 
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The  creation  of  two  incoming  particles  has  only  an  O(e)  success  probability,  which  can  be 
compensated  for  by  repeated  attempts.  (See  the  discussion  following  (f82ji.)  The  total  complexity 
of  preparing  two  particles  is  the  cost  of  simulating  the  time  evolution  given  in  (1751)  a  total  of  1/e 
times.  Thus,  by  the  results  of  113.31  the  complexity  is  (-^)1+°^.  Thus,  since  p  ~  a-1  for  fixed  e 
and  a  ~  e  for  fixed  p,  the  number  of  quantum  gates  Gexcite  needed  to  excite  the  two  initial  particles 
is 


|  £  3  ,  as  e  — >•  0  , 

\  p4+°(1)  ,  as  p  — >•  oo  . 


(74) 


Note  also  that  for  the  initial  wavepackets  to  be  well  separated,  L  must  be  larger  than  21.  Hence, 
in  the  high-momentum  limit  L  ~  p,  which  affects  the  complexity  of  other  steps  of  the  algorithm. 


Perturbative  Expansion 

The  resonant  excitation  can  be  analyzed  with  time-dependent  perturbation  theory.  Let 

1. 


-i  [  df  (H  +  XW(t)) 

Jo 


R  =  T  |  exp 

where  T{-}  denotes  the  time-ordered  product,  H  is  given  by 

W{t)=  f  dx  (7(i,x)Vv(x)  T  /*(f,x)V>jjQ(x))  , 


(75) 


(76) 


i  and  a  are  chosen  according  to  the  desired  type  of  particle,  and  /(t,x)  is  an  oscillatory  function 
whose  form  we  optimize  in  the  next  subsection.  The  end  product  of  the  excitation  process  is  i?|vac). 
One  can  expand  this  quantity  using  the  Dyson  series,  as  follows: 


R  =  1- 

where 

and  the  nth-order  term  in  A  is 


iX  [  dtiW}(ti)  +  (-iA)2  [  d  R  I"  dt2WI(t1)WI(t2) 
Jo  Jo  Jo 

W/(t)  =  eiHtW{t)e~im 


+ 


fT  ftn-1 

(-i\)n  dh...  dtnWI(t1)...WI(tn). 

Jo  Jo 

The  total  contribution  from  orders  A2  and  higher  is  bounded  by 

00  fT  ftn-1 

V(-iA)n  /  dh  . . .  /  dtnWj{U  ) . . .  W!{tn) 

Jo  Jo 


E 


A  nTr 


nl 


-w 


n= 2 

expfArrc]  —  1  —  A tw  , 


where 


w  =  max  ||W(i)||  . 

0<t<T 


(77) 

(78) 

(79) 

(80) 
(81) 

(82) 


From  the  above  analysis,  one  sees  that  the  Dyson  series  converges  rapidly.  The  single-particle 
excitation  amplitude  is  of  order  A,  and  the  dominant  error,  other  than  non-excitation,  is  the  two- 
particle  excitation  amplitude,  which  is  of  order  A2.  Setting  the  two-particle  excitation  probability 
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to  e,  one  obtains  a  single-particle  excitation  with  probability  p\  ~  yfe,  and  non-excitation  with 
probability  on  the  order  of  1  —  y/e.  In  a  standard  scattering  simulation,  one  wishes  to  prepare 
as  an  initial  state  single-particle  excitations  at  two  spatially  separated  locations.  The  fraction  of 
simulations  in  which  this  is  achieved  (rather  than  one  or  both  particles  failing  to  be  produced)  is 
thus  of  order  p\  ~  e.  One  can  detect  such  instances  and  compensate  by  repeating  the  simulation 
0(1 /p\)  times  and  postselecting  the  instances  in  which  both  particles  were  produced. 

Next,  we  consider  the  first-order  excitation  amplitude  in  more  detail.  Let  p)  be  any  state 
with  total  momentum  p  and  energy  E  above  the  vacuum  energy,  so  that  P\E,p)  =  p\E,  p)  and 
H\E,p)  =  E\E,  p),  where  P  is  the  total  momentum  operator.  (Here,  we  rely  on  the  fact  that 
[H,P]  =  0.)  Then,  to  first  order  in  A,  by  (|77l)  and  (f78]h 

(E,  p|i?|vac)  ~  —  iX  f  dt(E,p\Wi(t)\v&c)  (83) 

Jo 

=  —  iX  f  dte~'lEt(E,  p|VL(f)|vac) .  (84) 

Jo 

Recalling  that  the  momentum  operator  is  the  generator  of  spatial  translations,  one  has  ipit0l(x)  = 
e*PxV’i,a( 0)e“*Px.  Thus,  to  first  order  in  A, 


(E,  p|R|vac)  ~  —  iX  /  dt 

Jo 

(85) 

(Here  we  have  used  .P|vac}  =  0.)  Defining  f(t,  x)  =  0  for  t  (/  [0,r],  we  can  extend  the  time 
integration  to  infinity  and  express  (E,  p|R|vac)  in  terms  of  /,  the  Fourier  transform  of  /.  For  our 
choice  of  /,  given  in  the  next  subsection,  /  is  real,  and  therefore 


/ 


dxe-i(pt+px) 


/(L  x) (E,  p\ipi:a  (0) |vac)  +  /*(i,x)(£,p|^Ja(0)|vac 


(E,p\R\vac)  =  -iX  f(E,p)(E,p\^ita(0)\vac)  +  f(-E,-p)(E,p\tya(0)\vax;)  +  0( A2).  (86) 


Wavepacket  Shaping 


We  now  show  that  a  Gaussian  wavepacket  is  a  good  choice  for  /(t,x).  Specifically,  for  chosen 
constants  a,  /3  >  0,  let 


fit  x)  =  /  7/6X15  _  “  iEot  +  *Pox]  >  -r/2  <  t  <  t /2,-l/2  <  x  <  1/2,  ,g7. 

’ X  |  0 ,  otherwise .  ^ 


(For  convenience,  we  have  shifted  the  origin  of  the  coordinate  system.)  Here  rj  is  a  normalization 
factoid  with  mass  dimension  3/2.  With  this  choice  of  /, 


where 


f(E,  p)  =  r]qp,i(P  ~  Po)qa,r(E  -  E0) , 

rr/2 


Qp, 


rr/z 

■(d)  =  /  ■ 

J-r/2 


d  xeax-(px). 


(88) 


(89) 


3It  is  reasonable  to  choose  r/  so  that  fj  dtWi(t) |vac)  is  a  normalized  state.  In  the  ultrarelativistic  limit  this  implies 

that  rj  ~  (a2/?3 4  +  a4/32)1/4. 
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In  the  limit  r  — >•  oo,  the  function  qPyr{d)  converges  to  a  Gaussian  peak  of  width  ~  1  / p.  Since  E 
must  be  positive,  the  f(—E,  —  p)(E,  p|V^Q(0)lvac)  term  in  (l86|)  is  exponentially  small.  Hence,  one 
obtains 

(E,  p|R|vac)  ~  -iXf(E,  p )(E,  p| ^»,a(0)| vac).  (90) 

for  E  S>  1/r  and  A<  1.  By  (1901)  and  (1271).  one  sees  that  i?|vac)  is  a  antifermion  wavepacket  with 
momentum  centered  around  p.  To  create  a  fermion,  one  interchanges  ifj  and  ijfi  in  ()76l). 

Using  the  asymptotics  of  error  functions,  we  can  furthermore  bound  the  contributions  due  to  r 
being  finite.  One  finds  that 

| qPAd)  ~  Qp,oo(d)\  <  ^ e_(pr)2/4.  (91) 


4.4  Measuring  Number  Operators 

Recall  from  H3.ll  that  the  free  theory  (<?q  =  0)  is  exactly  solvable,  with  the  number  operators 
L-1aj (p)oj(p)  counting  fermions  of  species  j  in  momentum-mode  p  and  L~1bj(p)bj(p)  similarly 
counting  antifermions.  Thus,  as  one  possible  set  of  measurements  to  perform  on  the  final  state  of 
the  simulation,  we  propose,  as  in  [8j,  adiabatically  returning  to  the  free  theory  and  then  measuring 
number  operators  via  the  phase-estimation  algorithm.  We  analyze  this  measurement  procedure 
in  this  section.  An  alternative  set  of  measurements  that  is  more  suitable  when  bound  states  are 
present  is  analyzed  in  H4.5I 

The  adiabatic  return  to  the  free  theory  is  performed  in  the  presence  of  particle  wavepackets,  so 
the  state  being  adiabatically  transformed  is  not  an  energy  eigenstate.  Different  energy  eigenstates 
in  the  superposition  will  acquire  different  dynamical  phases  during  the  adiabatic  process  and  thus, 
in  physical  terms,  the  simulated  particles  will  propagate.  Such  propagation  is  undesirable  because 
we  do  not  want  any  scattering  to  occur  while  the  interaction  is  being  turned  off. 

Hence,  we  apply  the  same  technique  proposed  in  0  to  suppress  particle  propagation:  we  inter¬ 
leave  (simulated)  backwards  time  evolutions  governed  by  time-independent  Hamiltonians  into  the 
adiabatic  process.  By  an  analysis  similar  to  that  performed  in  [8],  one  finds  that,  to  ensure  that  a 
particle  propagates  no  further  than  a  distance  V,  it  suffices  to  use 

J  =  d{^)  <92> 

backwards  evolutions,  where  r  is  the  duration  of  the  original  adiabatic  process  and  p  is  the  mo¬ 
mentum  of  the  particle.  Further,  one  finds  that  the  total  probability  of  diabatically  exciting  one  or 
more  particles  i^ 

-Pdiabatic  =  o  ^ '  (93) 

Hence,  setting  V  to  a  constant  Pdiabatic  to  e,  one  obtains 

T  =  O  (^j  .  (94) 

4This  result  is  based  on  the  adiabatic  criterion  of  [35]  which  appears  to  be  applicable  [8]  to  our  Hamiltonian 
although  it  may  not  apply  to  all  Hamiltonians. 
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A  process  of  this  duration  can  be  implemented  with  (113.31) 


Gturn— off  —  0 


/  l2\  1+o(1)' 


- 


\ae  J 


(95) 


quantum  gates. 

The  phase-estimation  algorithm  [36]  enables  one  to  measure  in  the  eigenbasis  of  L~1Oj(p)aj(p), 
provided  one  can  efficiently  implement  e~lL  for  various  t  using  quantum  circuits.  By 

(1291)  and  (1301).  one  sees  that  the  problems  of  simulating  e~lL  aj(p)a-;(ph  and  its  antifermionic 
counterpart  are  largely  similar  to  the  problem  of  simulating  the  time  evolution  e~lHt,  which  was 
analyzed  in  detail  in  33.31  However,  these  number  operators  are  spatially  nonlocal,  which  means 
that  the  methods  of  33.31  do  not  perform  well  as  a  function  of  L.  Instead,  it  is  more  efficient  to  use 
recent  techniques  from  [37]. 

In  [37],  a  method  is  described  for  simulating  sparse  Hamiltonians  in  which  the  matrix  elements 
are  given  by  an  oracle.  As  discussed  on  pg.  2  of  [37] .  in  the  case  where  the  sparse  Hamiltonian 
consists  of  a  sum  of  d  terms  each  acting  on  0(1)  qubits,  the  number  of  oracle  queries  and  non-oracle- 
related  quantum  gates  both  scale  as  0{d).  A  number  operator  for  a  momentum  mode  consists  of 
0(L 2)  terms,  acting  between  all  pairs  of  spatial  lattice  sites.  Thus,  if  one  ignored  the  fermionic 
statistics,  the  number  of  non-oracle-related  gates  needed  to  simulate  the  time-evolution  induced 
by  a  number  operator  would  be  0(L2n)  =  0(L3).  The  number  of  gates  needed  to  implement  one 
oracle  query  to  the  sparse  matrix  defined  by  the  number  operator  would  be  O(n),  and  number  of 
quantum  gates  needed  to  implement  all  of  the  oracle  queries  would  be  0(L3).  Using  the  Bravyi- 
Kitaev  encoding  for  fermionic  statistics  adds  a  logarithmic  factor  to  the  complexity.  Measuring  all 
2NL  of  the  number  operators  thus  has  total  complexity  0(L 4)  =  0(L4/a4). 

4.5  Measuring  Local  Charge 

In  previous  work  [8],  we  proposed  measuring  local  energy  observables  as  an  alternative  to  returning 
to  the  free  theory  and  measuring  number  operators.  This  procedure  has  the  advantage  that  it  can 
detect  bound  states.  It  has  the  disadvantage  that  the  local  energy  observables  have  ultraviolet- 
divergent  vacuum  fluctuations  that  represent  a  noise  background  above  which  particle  excitations 
must  be  discerned.  In  this  paper,  we  instead  measure  simpler  local  observables,  namely  charges, 
whose  vacuum  fluctuations  are  less  difficult  to  control.  These  observables  can  thus  detect  charged 
bound  states,  although  they  are  blind  to  neutral  ones. 

From  the  equation  of  motion  of  the  massive  Gross-Neveu  model,  one  finds  that  for  each  j  € 
{1,  2, . . .  ,N}  the  quantity 

J f  (®)  =  ij>j  (x)7  (x)  (96) 

obeys 

=  0.  (97) 

Hence, 

Qj  =  Jj  =  & (x)7 Vj  (x)  (98) 

X  X 

is  a  conserved  charge.  Note  that,  for  any  b,  c  €  M,  Qj  =  bQj  +  c  is  also  conserved.  We  can  calibrate 

the  charge  observable  by  demanding  that  the  vacuum  have  zero  charge  and  that  particle  creation 
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change  the  charge  by  ±1.  One  satisfies  these  criteria  with  the  following  definition: 

Qj  =  «^j(x)70Vh‘(x)  -  Lt .  (99) 

By  (1271) ,  (1281) ,  and  (1361) ,  one  finds  that 

Qj  =  j  (aJ(p)ai(P)  “  b](p)bj(p))  ■  (10°) 

per 

For  any  envelope  function  /  :  0  — >•  [0, 1],  one  can  similarly  define 

Qjf)  =  /(x)  (a^i(x)7Vi(x)  -  1)  •  (101) 

( f) 

If  /  has  support  only  in  some  region  EcH,  then  Q-  can  be  thought  of  as  an  observable  for  the 
charge  in  that  region. 

The  most  obvious  choice  of  /  is  a  square  function  that  is  equal  to  one  inside  R  and  zero 
elsewhere.  However,  a  better  signal-to-noise  ratio  can  be  obtained  by  choosing  /  to  decay  from  one 
to  zero  more  smoothly  at  the  boundary  of  R.  Specifically,  calculations  (in  Appendix  |A])  show  that, 
when  /  is  chosen  to  be  a  Gaussian  of  width  R,  the  variance  of  the  observable  Q-  in  the  vacuum 
state  is  0(l/mR),  independent  of  the  lattice  spacing  a.  Hence  the  noise  background  above  which 
particle  excitations  are  to  be  detected  is  nondivergent  in  a  and  can  be  brought  to  an  arbitrarily 
low  level  at  the  cost  of  increasing  the  detector  size.  In  practice,  one  will  use  a  truncated  Gaussian, 
replacing  the  exponentially  small  tails  with  zero  at  distances  greater  than  some  constant  multiple 
of  R.  This  modified  /  then  has  support  on  a  region  of  size  O(R),  but  the  corresponding  operator 
is  exponentially  close  to  the  Gaussian  case  treated  by  our  analysis. 

Qj  has  eigenvalues  with  0(1)  separations.  Thus,  measuring  Qj  1  by  phase  estimation  entails 

simulating  the  unitary  transformation  exp  for  t  of  order  one.  Because  Q^'1  is  the  sum  of 

local  terms,  these  unitary  transformations  can  be  implemented  by  techniques  similar  to  those  in 
03. 31  with  complexity  C^gT1^1^-^1)). 

5  Some  Field-Theoretical  Aspects 

This  section  describes  some  quantum  field-theoretical  calculations:  analysis  of  the  effect  of  discretiz¬ 
ing  the  spatial  dimension  of  the  massive  Gross-Neveu  model,  and  the  perturbative  renormalization 
of  the  mass  in  the  discretized  theory. 

In  our  complexity  analysis  (112. 31).  our  criterion  for  choosing  the  lattice  spacing  a  is  that  the 
scattering  cross  sections  for  processes  at  a  momentum  scale  p  in  the  discretized  theory  should  differ 
from  their  continuum  values  by  at  most  a  factor  of  (1  +  e).  The  results  of  1)5.11  show  that  one 
can  satisfy  this  criterion  by  choosing  a  ~  e/p.  This  choice  then  affects  the  overall  scaling  of  the 
algorithm  in  the  large- momentum  and  high-precision  limits.  As  one  would  expect,  higher  energies 
and  greater  precision  require  a  smaller  lattice  spacing  and  thus  a  larger  number  of  lattice  sites 
(for  fixed  L).  Consequently,  the  number  of  quantum  gates  needed  to  simulate  time  evolutions  via 
Suzuki- Trotter  formulae  is  larger. 

In  H5.21  we  perturbatively  calculate  the  relationship  between  the  bare  mass  uiq,  which  is  a 
parameter  in  the  lattice  Hamiltonian  (see  (1221)  and  (1231)1.  and  the  physical  mass  m  of  the  particles 
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in  the  theory.  We  need  to  know  the  behavior  of  m  in  order  to  design  and  analyze  the  procedure 
for  preparing  the  interacting  vacuum  (H4.2J).  In  particular,  a  suitable  adiabatic  path  must  maintain 
a  non-zero  mass,  the  magnitude  of  which  affects  the  algorithmic  complexity,  as  indicated  by  the 
adiabatic  theorem. 

5.1  Effects  of  Non-zero  Lattice  Spacing 

The  effects  of  a  non-zero  lattice  spacing  can  be  analyzed  via  effective  field  theory.  The  discretized 
Lagrangian  can  be  thought  of  as  the  leading  contribution  to  an  effective  field  theory,  neglected 
terms  of  which  correspond  to  discretization  errors.  Hence,  the  scaling  of  the  error  with  the  lattice 
spacing  is  given  by  the  scaling  of  the  coefficients  of  those  terms. 

The  symmetries  of  the  continuum  theory  restrict  the  possible  operators  in  the  effective  held 
theory.  Consider  the  discrete  transformations  parity  (denoted  P ),  time  reversal  (T),  and  charge 
conjugation  (C).  Parity  changes  the  handedness  of  space  and  hence  reverses  the  momentum.  Thus, 


Pa(p)P  =  a(— p) , 

Pb(p)P  =  — 6(— p) . 

(102) 

IJsine  (HI)  and  (|102|).  we  then  obtain 

P^(t,x)P  =  7V(L-x) 

P#£,x)P  =  #£,-x)  7°. 

(103) 

Likewise, 

Ta(p)T  =  a(-p) , 

It  turns  out  that  time  reversal  needs  to  be  an 

T6(p)T  =  -6(-p). 

antilinear  operator.  Then 

(104) 

Tip(t,  x)T  =  7  V(-L  x) ,  Tip(t,  x)T  =  x)7X . 

Finally,  charge  conjugation  interchanges  particles  and  antiparticles.  Thus, 

(105) 

Ca{p)C  =  6(p) 

,  Cb(p)C  =  a{ p) , 

(106) 

and 

Cip(t,  x)C  =  ip*(t,  x) , 

Ci/j(t,x)C  =  ipT(t,x.) 7°  . 

(107) 

One  can  verify  that  the  Lagrangian  ([2])  is  invariant  under  each  of  the  transformations  P ,  T  and  C. 
Now  consider  the  operator  'ift'Wl'ift,  where  M  is  Hermitian.  Invariance  under  P,  T  and  C  requires 

M  =  7°M7°,  (108) 

M  =  -71M*71 ,  (109) 

M  =  -Mt.  (110) 

These  conditions  imply  that 

M  =  C7°,  c€l.  (Ill) 

Likewise,  for  where  M  is  Hermitian,  P,  T  and  C  invariance  requires 

M  =  (-l^My0, 

M  =  -(-i)VmV, 

M  =  Mt. 


22 


(112) 

(113) 

(114) 


These  conditions  imply  that,  for  /i  =  0, 

M  =  d  =  c(70)2,  cel,  (115) 

while,  for  fi=  1, 

M  =  cry5  =  — C7°71 ,  c  G  M  .  (116) 

Thus,  the  only  P-,  T-  and  C-invariant  bilinears  of  Dirac  fields  are  and  iif) {p,  =  0  or  1). 

Now  consider  four-fermion  operators,  namely,  products  of  two  bilinears.  The  set  {1  ,cr1}  forms 
a  complete  basis,  elements  of  which  satisfy  the  identity 

1  3 

SapfiyS  =  ^  ^  aaSajp)  ■  (H7) 

i=l 

For  7°  =  a2,  71  =  — iu1,  7 5  =  u3,  this  is  equivalent  to 

fiap3 7<5  =  2  T  (7^) a<5 (7 ^t) 7 /3  T  (7  )yp)  ■  (H8) 

Equation  (11181)  can  be  used  to  obtain  Fierz  identities.  For  example, 

=  {'^i)a(ipj)p{'ipj)'Y('ipi)sSap5^s 

=  ,  (119) 

where  the  minus  sign  comes  from  fermion  anticommutation.  Thus,  any  operator  of  the  form 
can  be  rewritten  as  a  sum  of  operators  of  the  form  7/i , F 1  f/7 ipj T 2 ipj ,  with  Fi2  G 

{1,7^,75}. 

If  Ti  7^  r 2 ,  then  will  violate  at  least  one  of  the  discrete  symmetries.  Furthermore, 

the  0(N )  symmetry^  associated  with  the  N  fermion  species  restricts  the  allowed  form  of  operators 
to  functions  of  YliLi  For  *  /  h  '&i'y5‘lPi‘ipjl5'll,j  is  ruled  out  by  invariance  under  P  (or  C) 

of  any  single  field  and  thus  (YliLi^Pn5^)  is  ruled  out.  Likewise,  (*  /  j)  and 

consequently  (  YliLi  are  ruled  out. 

We  conclude  that  the  only  four-fermion  operator  (without  derivatives)  in  the  effective  field 
theory  is  (^1 

Each  extra  derivative  or  factor  of  t/jTtl’  in  an  operator  will  increase  its  mass  dimension  by  one; 
correspondingly,  it  will  be  suppressed  by  an  extra  power  of  a.  We  therefore  conclude  that  no 
new  unsuppressed  operators  are  induced  in  the  effective  field  theory.  The  spatial  derivative  in  the 
continuum  theory  is  approximated  by  a  difference  operator,  with  an  error  of  order  a,  and  the  Wilson 
term  is  also  formally  of  order  a.  Spatial  discretization  errors  are  hence  of  order  a. 

5.2  Mass  Renormalization 

In  this  subsection,  we  calculate  the  renormalized  (or  physical)  mass  of  the  discretized  theory,  using 
second-order  perturbation  theory.  A  convenient  way  to  obtain  a  suitable  expression  is  to  use  a 
partially  renormalized  form  of  perturbation  theory  (as  was  done  in  [8] ) ,  in  which  one  uses  the  bare 
coupling  but  the  renormalized  mass. 

5In  fact,  the  massive  Gross-Neveu  model  has  an  0(2N )  symmetry. 
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To  perform  perturbative  calculations,  we  need  the  Feynman  rules  for  the  discretized  theory. 
The  propagator  is 

 +  fhjp) 


p2  —  m{jp ) 


2  5 


where 


/  o  1  ,  ,  .  2r  .  2  ( apl\ 

P  =  P  ,  ~  sm(ap  )  ,  m(p)  =  m  H - sm  - 

V  a  )  a  V  2  / 


(120) 


(121) 


For  convenience,  we  use  the  standard  technique  of  introducing  an  auxiliary  held  a  and  rewrite  the 
Lagrangian  as 

C  =  Co  +  Ca  ,  (122) 


where  Cq  is  the  discretized  free  Lagrangian  and 


£CT  =  -^(T2  -  ga'ipj'i/jj 


The  corresponding  Feynman  rules  are 


=  —  i . 


=  -ig- 


(123) 


(124) 


At  one-loop  order, 


—  iM  ( p )  = 


+ 


where  the  second  diagram  is  the  counterterm. 
The  first  diagram  gives 


=  -So 


dk°  r'a  dk 1  +  fn(k) 


r°°  dk u  r/a 

J—o o  2vr  J—ir/a 

/IT 

dk 1 

-7 r 


-tt/o  2vr 

ma  +  2r  sin2  (^-) 

-ysin2  L1  +  (ma  +  2rsin2  (^-))2 


(125) 


(126) 

(127) 


The  term  in  (|127D  proportional  to  r  scales  as  1/a  and  gives  the  mass  correction  to  the  doubler 
(spurious  fermion).  The  term  proportional  to  m  gives  the  following: 


9omi  /  \  . 

mo  =  m - - —  iog(maj  +  •  •  •  . 

Z7T 

At  two-loop  order,  the  1PI  amplitude  has  the  additional  contributions 

,0\ 


(128) 


+ 


+ 


+ 


The  renormalization  condition  satished  at  hrst  order  implies  that  the  first  two  diagrams  cancel. 
The  last  two  diagrams  give 


■  4 

m 

167T3 


-I 

a 


(129) 
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and 


,0\ 


igtN 

167T3 


+  i/f 

a 


(130) 


where  I |a^ ,  if^ ,  if^  and  if'*  are  given  in  Appendix  [Bj  Numerical  evaluation  of  these  integrals 
reveals  the  forms 


if ma  +  •  •  •  ,  (131) 

if^  =  cfV  log2  (ma)  -  cf2)  log  (ma)  +  cf^  H - ,  (132) 

(i) 

with  c)  >  0.  We  thus  obtain 

4  (1) 

m  =  (N cf1'>  —  cf1'1 )  log2 (m^  a)  +  ■  ■  ■  ,  (133) 

lovr-5  ' 


where  m ^  denotes  the  physical  mass  at  one-loop  order. 
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A  Variance  of  Local  Charge 

Consider  Qy-  in  the  continuum  limit: 

Q{P  =  j  dxjf(x)/(x).  (134) 

We  wish  to  compute  its  variance,  which  is  given  by 

({QP  ~  (Qp))2)  =  J  dxdyf(x)f(y)({Jj°\x.)Jj°\y))  —  (jj0)(x))(jj0)(y)))  (135) 

/rl2k  ~  ~ 

—  \f(k)\2Gc(k),  (136) 


where  Gc  is  the  connected  Green’s  function.  By  standard  quantum  field-theoretical  methods,  we 
obtain 


Gc{k°,kl) 


f1  f  d2pE  m2  -  x(l  -  x)((k0)2  +  (k1)2) 

1 J0  J  (27t)2  \p2E  +  m2  —  x(\  —  x)((fc0)2  —  (A;1)2)]2  ’ 
i  f1  m2  —  x{l  —  x)((k0)2  +  (k1)2) 

27 r  J0  X  m?  —  x(l  —  x)((k0)2  —  (k1)2)  ' 


Substituting  Gc  into  (I136D  and  using  an  ultraviolet  regulator,  we  obtain 


{QP  -  (Q 


(/) 


))2  = 


(2vr) 


dktlfik1)]2  J  dx 


(k1) 


^)2  +  ^ 


(137) 

(138) 


(139) 


For  the  square  window  function  /(x)  =  9(R/2  —  |x|),  the  charge  fluctuation  diverges,  but  only 
logarithmically,  with 


{[QP  ~  { QP ))2)  <  ^(log(27r)  -  1  -  log(ma))  +  0((ma)2) . 
For  /(x)  =  exp(— x2/R2), 


{(qt  -  <<?r»2} 


V^vr3/2  1 

32  mR  +  "  '  ' 


(140) 


(141) 
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B  Integrals  for  Mass  Renormalization 


For  i  =  1,2, 


t(°) 

1i 


Ab) 

1,  =  r 


1  j  j  j  <J(x  +  y  +  «-  1) 
ax  ay  az  — -=  . 

o  y/xy  +  xz  +  yz 


where 


JV- 


JV- 


N- 


N- 


N< 


N. 


N. 


N. 


l 


,  ,  ,  <5(x  +  y  +  z-  1) 

ax  ay  az  — -= — ^=^=-  , 

o  y/xy  +  xz  +  yz  J_ 


/7T  /*7 T  / 

jkLd«  ( 

/7T  /*7T 

dk  /  hy 

-7T  J  — 7T 


iV> 


(al) 


2(xy  +  xz  +  yz)2D 


+ 


N„ 


(61) 


2(xy  +  xz  +  yz)2D 


D  =  x 


(al) 


(a2) 


sin2  (A;)  +  (  ma  +  2r  sin"  I  — 


A;\  \  2i 


+2 


+  y 

q  —  k\  \  2i 


sin2(y  —  A;)  +  |ma  +  2r  sin2  ^ 

cyz  — 

2  lk~Q 


xyz 


xy  +  xz  +  yz 


(ma) 


9  0  0  0  0  0 
x  y  —  xy  +  x  z  +  4 xyz  —  y  z  +  xz  +  yz 


4/3r  sin"  |  -  )  sin' 


+  /3  sin(Ai)  sin  (A:  —  y)  +  4/2r2  sin" 


]y(a2)  \ 

F>2  J 

,  (142) 

A7.(62)  \ 

Z?2  j 

,  (143) 

D)1 

(144) 

(145) 

-■  (i: 

(146) 

+4/ir2  sin2  (yy  )  sin2  (|)  ~  ^2  sin^^  +  ^  sin^  “  ?)  sin(^) 
+2(mfl)r|/2/3  sin2  [^j  +  /i/3  sin2  (yy)  +  f\h  sm2  ( | ) )  +  (ma)2/i/2/3  , 


/i  = 
(61)  _ 


xy  +  xz  +  2yz 


2xy  +  xz  +  yz 

J2  |  |  5  j3 


(62)  _ 


xy  +  xz  +  yz  "  “  xy  +  xz  +  yz 

x2y  —  xy2  +  x2z  +  xyz  —  y2z  +  xz2  +  yz2 
—  (xy  +  xz  +  yz)(x cos(fc)  —  y  cos(y)  +  z cos(/c  —  y)) , 
k 

2)  V  2 

2  I k  —  q\  .  ,  ,  ~  .  9 (k 


xy  +  xz 
xy  +  xz  +  yz 


q2  ■  2  (  ^  \  ■  2  (  k  Q\  ■  2  { Q 
8 r  sin  —  sin  '  1  1 


(al)  _ 


(a2)  _ 


—2  sin  (A:)  sin2 
4z(xy  +  xz  +  yz) , 

t  ^2  (^)2 

—  (ma)  t - ttt 

(xy  +  xz  +  yz)-3 

2xy  +  xz  +  yz 


sin2  +  2sin(A;)  sin(A;  —  q)  sin2 
sin(y)  +  2  sin2  ^  sin(A:  —  q)  sin(y) , 

xyz2  [ma  +  2 r  sin2  [j-^j  ) 


—  (ma)- 


(xy  +  xz  +  yz)2 


(147) 

(148) 

(149) 

(150) 


+ 


xy  +  xz  +  yz 


^  |ma  +  2r  sin2  [— ^  [ma  +  2 r  sin2  j  —  sin(A:)  sin(y)^ 


(61)  _ 


+2r  sin2  (^)  (ma  +  2r(  sin2  f +  sin2  ( ? 


.  2  fk-q 


(62)  _ 


2 z(xy  +  xz  +  yz)  sin 

^  ^4r2  sin2 


.  9  (k  —  q 
2  sin2  1  1 


2 

.  2  ■  n\  ■ 

sm  (  —  )  —  sm[k)  sin 


(151) 

(152) 
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